Exploration into trying to detect climate regions using the rainfall data.

Initialization

library(knitr)
source("./Scripts/Settings.R")
source("./Scripts/AusWeather.R")
source("./Scripts/AusClimate.R")
source("./Scripts/AusMaps.R")
source("./Scripts/Utility.R")

Load rainfall valuations

regional_valuations <- function() {
    message("Evaluating Baseline Stations")
    baseline_stations <- aus_baseline_stations(aus_rainfall_valuation_info)
    message("Loading Observations")
    observations <- aus_clean_observations(baseline_stations$Station, aus_rainfall_clean_info)
    message("Evaluating Baselines")
    baselines <- aus_valuation_baselines(observations, aus_rainfall_valuation_info)
    message("Evaluating Valuations")
    valuations <- aus_valuations(observations, aus_rainfall_valuation_info, baselines = baselines)
    message("Evaluation Complete")
    list(
        stations = baseline_stations,
        observations = observations,
        baselines = baselines,
        valuations = valuations
    )
}
memorised_regional_valuations <- aus_memorise(regional_valuations, aus_rainfall_observations_info,
                                              paste0("RainRegionsValuations v", aus_rainfall_version))
stations <- memorised_regional_valuations$stations
observations <- memorised_regional_valuations$observations
baselines <- memorised_regional_valuations$baselines
valuations <- memorised_regional_valuations$valuations

Regional Exploration

Baseline

Plot a map of australia wide rainfall baselines for each month

baselines %>%
    inner_join(stations, by = c("Zone", "Station")) %>%
    ggplot() %>%
    aus_dark_map() + 
    geom_point(aes(x = Longitude, y = Latitude, color = log(Baseline.Mean)), size = 0.1) +
    scale_colour_distiller(type = "div", palette="RdBu", direction = 1) +
    guides(colour=FALSE) +
    facet_wrap(~ Minor_Period)

NA

Try using kmeans to try to detect regions with similar variation.

Anaomly Map

Plot a map of australia wide rainfall anamolies for a series of years.

valuations %>%
    filter(Major_Period >= 2010, Major_Period <= 2015, Minor_Period == 1) %>%
    inner_join(stations, by = c("Zone", "Station")) %>%
    ggplot() %>%
    aus_dark_map() + 
    geom_point(aes(x = Longitude, y = Latitude, color = Rainfall.Anamoly), size = 0.1) +
    scale_colour_distiller(type = "div", palette="RdBu", direction = 1) +
    guides(colour=FALSE) +
    facet_wrap(~ Major_Period)

It looks like anamolies occur in distinct regions. And its quite easy to see.

Try using kmeans to try to detect regions that vary together. Also could use correlation. The corclust method in the klaR package can be used to detect clusters.

Try Regional Clustering

Try finding clusters using baseline rainfalls. Spread the monthly rainfalls into columns

station_rainfall_means <- baselines %>% 
    select(Zone, Station, Rainfall.Mean = Minor_Period, Baseline.Mean) %>% 
    spread(key = Rainfall.Mean, value = Baseline.Mean, sep = ".") %>%
    ungroup()
station_rainfall_sds <- baselines %>% 
    select(Zone, Station, Rainfall.SD = Minor_Period, Baseline.SD) %>% 
    spread(key = Rainfall.SD, value = Baseline.SD, sep = ".") %>%
    ungroup()
station_rainfall <- station_rainfall_means %>% inner_join(station_rainfall_sds, by = c("Zone", "Station"))
station_rainfall

Group stations by means

station_centers_mean_groups <- function(centers) {
    station_kmeans <- station_rainfall_means %>% select(-Zone, -Station) %>% kmeans(centers)
    station_rainfall_means$Centers <- centers
    station_rainfall_means$Group <- station_kmeans$cluster
    station_rainfall_means
}
centers <- 4:20
station_rainfall_mean_by_centers <- centers %>% map_dfr(station_centers_mean_groups)
for (center in centers) {
    p <- station_rainfall_mean_by_centers %>%
        inner_join(stations, by = c("Zone", "Station")) %>%
        filter(Centers == center) %>%
        ggplot() %>%
        aus_dark_map() + 
        geom_point(aes(x = Longitude, y = Latitude, color = Group)) +
        scale_colour_distiller(type = "div", palette="RdBu") +
        guides(colour=FALSE) +
        ggtitle(paste0("Centers = ", center))
    print(p)
}

Group stations by sds

station_centers_sd_groups <- function(centers) {
    station_kmeans <- station_rainfall_sds %>% select(-Zone, -Station) %>% kmeans(centers)
    station_rainfall_sds$Centers <- centers
    station_rainfall_sds$Group <- station_kmeans$cluster
    station_rainfall_sds
}
centers <- 4:20
station_rainfall_sd_by_centers <- centers %>% map_dfr(station_centers_sd_groups)
did not converge in 10 iterations
for (center in centers) {
    p <- station_rainfall_sd_by_centers %>%
        inner_join(stations, by = c("Zone", "Station")) %>%
        filter(Centers == center) %>%
        ggplot() %>%
        aus_dark_map() + 
        geom_point(aes(x = Longitude, y = Latitude, color = Group)) +
        scale_colour_distiller(type = "div", palette="RdBu") +
        guides(colour=FALSE)
    print(p)
}

Regional Grouping Conclusion

SDs alone it not so distinct and the algor has trouble. Use just the means so their is no issue in getting values on the same scale.

Finds group with similar rainfall distributions. Picks out the arid regions, coast regions and so one.

However tends to do things like group areas on the east and west coast of Aus which I know tend to vary seperately.

10 Centers looks about nice. It picks out the areas that are in drought right now.

Save as station_regional_clusters

station_regional_clusters <- station_rainfall_means %>% 
    select(-Zone, -Station) %>% 
    kmeans(10)
station_regions <- station_rainfall_means
station_regions$RegionID <- station_regional_clusters$cluster
station_regions <- station_regions %>% select(Zone, Station, RegionID, everything())
stations %>%
    inner_join(station_regions, by = c("Zone", "Station")) %>%
    ggplot() %>%
    aus_dark_map() + 
    geom_point(aes(x = Longitude, y = Latitude, color = RegionID))  +
    scale_colour_distiller(type = "div", palette="RdBu")

Regional Group Anamoly

Generate anamoly grouped by regions

Combined plot

plot(valuations %>%
    inner_join(stations, by = c("Zone", "Station")) %>%
    inner_join(station_regions, by = c("Zone", "Station")) %>%
    filter(Major_Period > 1900) %>%
    group_by(RegionID, Median_Date) %>%
    summarise(Rainfall.Anamoly = mean(Rainfall.Anamoly)) %>%
    ggplot(aes(x = Median_Date, y = Rainfall.Anamoly)) +
    geom_line(color="lightgrey") +
    geom_smooth(method = "loess", formula = y ~ x, span = .1) +
    facet_wrap(~ RegionID) + ylim(-1, 1))

Paged Plot

for (regionID in 1:10) {
    
    gridExtra::grid.arrange(
        stations %>%
            inner_join(station_regions, by = c("Zone", "Station")) %>%
            filter(RegionID == regionID) %>%
            ggplot() %>%
            aus_dark_map() + 
            geom_point(aes(x = Longitude, y = Latitude)) +
            guides(colour=FALSE),
        
        valuations %>%
            inner_join(stations, by = c("Zone", "Station")) %>%
            inner_join(station_regions, by = c("Zone", "Station")) %>%
            filter(RegionID == regionID) %>%
            filter(Major_Period > 1900) %>%
            group_by(Median_Date) %>%
            summarise(Rainfall.Anamoly = mean(Rainfall.Anamoly)) %>%
            ggplot(aes(x = Median_Date, y = Rainfall.Anamoly)) +
            geom_line(color="grey") +
            geom_smooth(method = "loess", formula = y ~ x, span = .2) +
             ylim(-1, 1),
        
        ncol = 2
    )
}

Climate Grouping Conclusion

It sorta works. You can see the major drought in eastern australia. However some areas are associated that probably shouldn’t be.

Try Correlative Clustering

Cluster stations that tend to vary together.

Try kmeans to start cos its easy!

Try spreading all periods during the baseline period into one set and use kmeans to group stations.

Group stations by anamoly

station_rainfall_anamoly_clusters <- function(centers) {
    station_kmeans <- station_rainfall_anamoly %>% select(-Zone, -Station) %>% kmeans(centers)
    r <- station_rainfall_anamoly
    r$Centers <- centers
    r$RegionID <- station_kmeans$cluster
    r <- r %>% select(Zone, Station, Centers, RegionID)
    r
}
anaomly_centers <- 4:30
station_rainfall_anamoly_centers <- anaomly_centers %>% map_dfr(station_rainfall_anamoly_clusters)
for (center in anaomly_centers) {
    plot(station_rainfall_anamoly_centers %>%
        inner_join(stations, by = c("Zone", "Station")) %>%
        filter(Centers == center) %>%
        ggplot() %>%
        aus_dark_map() + 
        geom_point(aes(x = Longitude, y = Latitude, color = RegionID)) +
        scale_colour_distiller(type = "div", palette="RdBu") +
        guides(colour=FALSE) +
        ggtitle(paste0("Centers = ", center)))
}

Paged Plot

selected_anaomly_centers <- 15
selected_station_rainfall_anamoly <- station_rainfall_anamoly_centers %>% filter(Centers == selected_anaomly_centers)
selected_anaomly_start <- 1950
selected_anaomly_end <- 2018
for (regionID in 1:selected_anaomly_centers) {
    
    gridExtra::grid.arrange(
        stations %>%
            inner_join(selected_station_rainfall_anamoly, by = c("Zone", "Station")) %>%
            filter(RegionID == regionID) %>%
            ggplot() %>%
            aus_dark_map() + 
            geom_point(aes(x = Longitude, y = Latitude)) +
            guides(colour=FALSE),
        
        valuations %>%
            inner_join(stations, by = c("Zone", "Station")) %>%
            inner_join(selected_station_rainfall_anamoly, by = c("Zone", "Station")) %>%
            filter(RegionID == regionID) %>%
            filter(Major_Period >= selected_anaomly_start, Major_Period <= selected_anaomly_end) %>%
            group_by(Median_Date) %>%
            summarise(Rainfall.Anamoly = mean(Rainfall.Anamoly)) %>%
            ggplot(aes(x = Median_Date, y = Rainfall.Anamoly)) +
            geom_line(color="grey") +
            geom_smooth(method = "loess", formula = y ~ x, span = .1),
        
        ncol = 2
    )
}

LS0tDQp0aXRsZTogIkF1c3RyYWxpYW4gQ2xpbWF0ZSAtIFJhaW5mYWxsIFJlZ2lvbnMiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpFeHBsb3JhdGlvbiBpbnRvIHRyeWluZyB0byBkZXRlY3QgY2xpbWF0ZSByZWdpb25zIHVzaW5nIHRoZSByYWluZmFsbCBkYXRhLg0KDQojIEluaXRpYWxpemF0aW9uDQoNCmBgYHtyIG1lc3NhZ2U9RkFMU0V9DQoNCmxpYnJhcnkoa25pdHIpDQoNCnNvdXJjZSgiLi9TY3JpcHRzL1NldHRpbmdzLlIiKQ0Kc291cmNlKCIuL1NjcmlwdHMvQXVzV2VhdGhlci5SIikNCnNvdXJjZSgiLi9TY3JpcHRzL0F1c0NsaW1hdGUuUiIpDQpzb3VyY2UoIi4vU2NyaXB0cy9BdXNNYXBzLlIiKQ0Kc291cmNlKCIuL1NjcmlwdHMvVXRpbGl0eS5SIikNCg0KYGBgDQoNCiMjIExvYWQgcmFpbmZhbGwgdmFsdWF0aW9ucw0KDQpgYGB7cn0NCg0KcmVnaW9uYWxfdmFsdWF0aW9ucyA8LSBmdW5jdGlvbigpIHsNCiAgICBtZXNzYWdlKCJFdmFsdWF0aW5nIEJhc2VsaW5lIFN0YXRpb25zIikNCiAgICBiYXNlbGluZV9zdGF0aW9ucyA8LSBhdXNfYmFzZWxpbmVfc3RhdGlvbnMoYXVzX3JhaW5mYWxsX3ZhbHVhdGlvbl9pbmZvKQ0KICAgIG1lc3NhZ2UoIkxvYWRpbmcgT2JzZXJ2YXRpb25zIikNCiAgICBvYnNlcnZhdGlvbnMgPC0gYXVzX2NsZWFuX29ic2VydmF0aW9ucyhiYXNlbGluZV9zdGF0aW9ucyRTdGF0aW9uLCBhdXNfcmFpbmZhbGxfY2xlYW5faW5mbykNCiAgICBtZXNzYWdlKCJFdmFsdWF0aW5nIEJhc2VsaW5lcyIpDQogICAgYmFzZWxpbmVzIDwtIGF1c192YWx1YXRpb25fYmFzZWxpbmVzKG9ic2VydmF0aW9ucywgYXVzX3JhaW5mYWxsX3ZhbHVhdGlvbl9pbmZvKQ0KICAgIG1lc3NhZ2UoIkV2YWx1YXRpbmcgVmFsdWF0aW9ucyIpDQogICAgdmFsdWF0aW9ucyA8LSBhdXNfdmFsdWF0aW9ucyhvYnNlcnZhdGlvbnMsIGF1c19yYWluZmFsbF92YWx1YXRpb25faW5mbywgYmFzZWxpbmVzID0gYmFzZWxpbmVzKQ0KICAgIG1lc3NhZ2UoIkV2YWx1YXRpb24gQ29tcGxldGUiKQ0KICAgIGxpc3QoDQogICAgICAgIHN0YXRpb25zID0gYmFzZWxpbmVfc3RhdGlvbnMsDQogICAgICAgIG9ic2VydmF0aW9ucyA9IG9ic2VydmF0aW9ucywNCiAgICAgICAgYmFzZWxpbmVzID0gYmFzZWxpbmVzLA0KICAgICAgICB2YWx1YXRpb25zID0gdmFsdWF0aW9ucw0KICAgICkNCn0NCg0KbWVtb3Jpc2VkX3JlZ2lvbmFsX3ZhbHVhdGlvbnMgPC0gYXVzX21lbW9yaXNlKHJlZ2lvbmFsX3ZhbHVhdGlvbnMsIGF1c19yYWluZmFsbF9vYnNlcnZhdGlvbnNfaW5mbywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwYXN0ZTAoIlJhaW5SZWdpb25zVmFsdWF0aW9ucyB2IiwgYXVzX3JhaW5mYWxsX3ZlcnNpb24pKQ0KDQpzdGF0aW9ucyA8LSBtZW1vcmlzZWRfcmVnaW9uYWxfdmFsdWF0aW9ucyRzdGF0aW9ucw0Kb2JzZXJ2YXRpb25zIDwtIG1lbW9yaXNlZF9yZWdpb25hbF92YWx1YXRpb25zJG9ic2VydmF0aW9ucw0KYmFzZWxpbmVzIDwtIG1lbW9yaXNlZF9yZWdpb25hbF92YWx1YXRpb25zJGJhc2VsaW5lcw0KdmFsdWF0aW9ucyA8LSBtZW1vcmlzZWRfcmVnaW9uYWxfdmFsdWF0aW9ucyR2YWx1YXRpb25zDQoNCmBgYA0KDQojIFJlZ2lvbmFsIEV4cGxvcmF0aW9uDQoNCiMjIEJhc2VsaW5lDQoNClBsb3QgYSBtYXAgb2YgYXVzdHJhbGlhIHdpZGUgcmFpbmZhbGwgYmFzZWxpbmVzIGZvciBlYWNoIG1vbnRoDQoNCmBgYHtyfQ0KDQpiYXNlbGluZXMgJT4lDQogICAgaW5uZXJfam9pbihzdGF0aW9ucywgYnkgPSBjKCJab25lIiwgIlN0YXRpb24iKSkgJT4lDQogICAgZ2dwbG90KCkgJT4lDQogICAgYXVzX2RhcmtfbWFwKCkgKyANCiAgICBnZW9tX3BvaW50KGFlcyh4ID0gTG9uZ2l0dWRlLCB5ID0gTGF0aXR1ZGUsIGNvbG9yID0gbG9nKEJhc2VsaW5lLk1lYW4pKSwgc2l6ZSA9IDAuMSkgKw0KICAgIHNjYWxlX2NvbG91cl9kaXN0aWxsZXIodHlwZSA9ICJkaXYiLCBwYWxldHRlPSJSZEJ1IiwgZGlyZWN0aW9uID0gMSkgKw0KICAgIGd1aWRlcyhjb2xvdXI9RkFMU0UpICsNCiAgICBmYWNldF93cmFwKH4gTWlub3JfUGVyaW9kKQ0KDQogICAgDQpgYGANCg0KVHJ5IHVzaW5nIGttZWFucyB0byB0cnkgdG8gZGV0ZWN0IHJlZ2lvbnMgd2l0aCBzaW1pbGFyIHZhcmlhdGlvbi4NCg0KDQojIyBBbmFvbWx5IE1hcA0KDQpQbG90IGEgbWFwIG9mIGF1c3RyYWxpYSB3aWRlIHJhaW5mYWxsIGFuYW1vbGllcyBmb3IgYSBzZXJpZXMgb2YgeWVhcnMuDQoNCmBgYHtyfQ0KDQp2YWx1YXRpb25zICU+JQ0KICAgIGZpbHRlcihNYWpvcl9QZXJpb2QgPj0gMjAxMCwgTWFqb3JfUGVyaW9kIDw9IDIwMTUsIE1pbm9yX1BlcmlvZCA9PSAxKSAlPiUNCiAgICBpbm5lcl9qb2luKHN0YXRpb25zLCBieSA9IGMoIlpvbmUiLCAiU3RhdGlvbiIpKSAlPiUNCiAgICBnZ3Bsb3QoKSAlPiUNCiAgICBhdXNfZGFya19tYXAoKSArIA0KICAgIGdlb21fcG9pbnQoYWVzKHggPSBMb25naXR1ZGUsIHkgPSBMYXRpdHVkZSwgY29sb3IgPSBSYWluZmFsbC5BbmFtb2x5KSwgc2l6ZSA9IDAuMSkgKw0KICAgIHNjYWxlX2NvbG91cl9kaXN0aWxsZXIodHlwZSA9ICJkaXYiLCBwYWxldHRlPSJSZEJ1IiwgZGlyZWN0aW9uID0gMSkgKw0KICAgIGd1aWRlcyhjb2xvdXI9RkFMU0UpICsNCiAgICBmYWNldF93cmFwKH4gTWFqb3JfUGVyaW9kKQ0KDQpgYGANCg0KSXQgbG9va3MgbGlrZSBhbmFtb2xpZXMgb2NjdXIgaW4gZGlzdGluY3QgcmVnaW9ucy4gQW5kIGl0cyBxdWl0ZSBlYXN5IHRvIHNlZS4NCg0KVHJ5IHVzaW5nIGttZWFucyB0byB0cnkgdG8gZGV0ZWN0IHJlZ2lvbnMgdGhhdCB2YXJ5IHRvZ2V0aGVyLiBBbHNvIGNvdWxkIHVzZSBjb3JyZWxhdGlvbi4gVGhlIGNvcmNsdXN0IG1ldGhvZCBpbiB0aGUga2xhUiBwYWNrYWdlIGNhbiBiZSB1c2VkIHRvIGRldGVjdCBjbHVzdGVycy4NCg0KIyBUcnkgUmVnaW9uYWwgQ2x1c3RlcmluZw0KDQpUcnkgZmluZGluZyBjbHVzdGVycyB1c2luZyBiYXNlbGluZSByYWluZmFsbHMuIFNwcmVhZCB0aGUgbW9udGhseSByYWluZmFsbHMgaW50byBjb2x1bW5zDQoNCg0KYGBge3J9DQoNCnN0YXRpb25fcmFpbmZhbGxfbWVhbnMgPC0gYmFzZWxpbmVzICU+JSANCiAgICBzZWxlY3QoWm9uZSwgU3RhdGlvbiwgUmFpbmZhbGwuTWVhbiA9IE1pbm9yX1BlcmlvZCwgQmFzZWxpbmUuTWVhbikgJT4lIA0KICAgIHNwcmVhZChrZXkgPSBSYWluZmFsbC5NZWFuLCB2YWx1ZSA9IEJhc2VsaW5lLk1lYW4sIHNlcCA9ICIuIikgJT4lDQogICAgdW5ncm91cCgpDQoNCnN0YXRpb25fcmFpbmZhbGxfc2RzIDwtIGJhc2VsaW5lcyAlPiUgDQogICAgc2VsZWN0KFpvbmUsIFN0YXRpb24sIFJhaW5mYWxsLlNEID0gTWlub3JfUGVyaW9kLCBCYXNlbGluZS5TRCkgJT4lIA0KICAgIHNwcmVhZChrZXkgPSBSYWluZmFsbC5TRCwgdmFsdWUgPSBCYXNlbGluZS5TRCwgc2VwID0gIi4iKSAlPiUNCiAgICB1bmdyb3VwKCkNCg0Kc3RhdGlvbl9yYWluZmFsbCA8LSBzdGF0aW9uX3JhaW5mYWxsX21lYW5zICU+JSBpbm5lcl9qb2luKHN0YXRpb25fcmFpbmZhbGxfc2RzLCBieSA9IGMoIlpvbmUiLCAiU3RhdGlvbiIpKQ0Kc3RhdGlvbl9yYWluZmFsbA0KDQoNCmBgYA0KDQojIyBHcm91cCBzdGF0aW9ucyBieSBtZWFucw0KDQpgYGB7cn0NCg0Kc3RhdGlvbl9jZW50ZXJzX21lYW5fZ3JvdXBzIDwtIGZ1bmN0aW9uKGNlbnRlcnMpIHsNCiAgICBzdGF0aW9uX2ttZWFucyA8LSBzdGF0aW9uX3JhaW5mYWxsX21lYW5zICU+JSBzZWxlY3QoLVpvbmUsIC1TdGF0aW9uKSAlPiUga21lYW5zKGNlbnRlcnMpDQogICAgc3RhdGlvbl9yYWluZmFsbF9tZWFucyRDZW50ZXJzIDwtIGNlbnRlcnMNCiAgICBzdGF0aW9uX3JhaW5mYWxsX21lYW5zJEdyb3VwIDwtIHN0YXRpb25fa21lYW5zJGNsdXN0ZXINCiAgICBzdGF0aW9uX3JhaW5mYWxsX21lYW5zDQp9DQoNCmNlbnRlcnMgPC0gNDoyMA0Kc3RhdGlvbl9yYWluZmFsbF9tZWFuX2J5X2NlbnRlcnMgPC0gY2VudGVycyAlPiUgbWFwX2RmcihzdGF0aW9uX2NlbnRlcnNfbWVhbl9ncm91cHMpDQoNCmZvciAoY2VudGVyIGluIGNlbnRlcnMpIHsNCiAgICBwIDwtIHN0YXRpb25fcmFpbmZhbGxfbWVhbl9ieV9jZW50ZXJzICU+JQ0KICAgICAgICBpbm5lcl9qb2luKHN0YXRpb25zLCBieSA9IGMoIlpvbmUiLCAiU3RhdGlvbiIpKSAlPiUNCiAgICAgICAgZmlsdGVyKENlbnRlcnMgPT0gY2VudGVyKSAlPiUNCiAgICAgICAgZ2dwbG90KCkgJT4lDQogICAgICAgIGF1c19kYXJrX21hcCgpICsgDQogICAgICAgIGdlb21fcG9pbnQoYWVzKHggPSBMb25naXR1ZGUsIHkgPSBMYXRpdHVkZSwgY29sb3IgPSBHcm91cCkpICsNCiAgICAgICAgc2NhbGVfY29sb3VyX2Rpc3RpbGxlcih0eXBlID0gImRpdiIsIHBhbGV0dGU9IlJkQnUiKSArDQogICAgICAgIGd1aWRlcyhjb2xvdXI9RkFMU0UpICsNCiAgICAgICAgZ2d0aXRsZShwYXN0ZTAoIkNlbnRlcnMgPSAiLCBjZW50ZXIpKQ0KICAgIHByaW50KHApDQp9DQoNCmBgYA0KDQoNCiMjIEdyb3VwIHN0YXRpb25zIGJ5IHNkcw0KDQpgYGB7cn0NCg0Kc3RhdGlvbl9jZW50ZXJzX3NkX2dyb3VwcyA8LSBmdW5jdGlvbihjZW50ZXJzKSB7DQogICAgc3RhdGlvbl9rbWVhbnMgPC0gc3RhdGlvbl9yYWluZmFsbF9zZHMgJT4lIHNlbGVjdCgtWm9uZSwgLVN0YXRpb24pICU+JSBrbWVhbnMoY2VudGVycykNCiAgICBzdGF0aW9uX3JhaW5mYWxsX3NkcyRDZW50ZXJzIDwtIGNlbnRlcnMNCiAgICBzdGF0aW9uX3JhaW5mYWxsX3NkcyRHcm91cCA8LSBzdGF0aW9uX2ttZWFucyRjbHVzdGVyDQogICAgc3RhdGlvbl9yYWluZmFsbF9zZHMNCn0NCg0KY2VudGVycyA8LSA0OjIwDQpzdGF0aW9uX3JhaW5mYWxsX3NkX2J5X2NlbnRlcnMgPC0gY2VudGVycyAlPiUgbWFwX2RmcihzdGF0aW9uX2NlbnRlcnNfc2RfZ3JvdXBzKQ0KDQpmb3IgKGNlbnRlciBpbiBjZW50ZXJzKSB7DQogICAgcCA8LSBzdGF0aW9uX3JhaW5mYWxsX3NkX2J5X2NlbnRlcnMgJT4lDQogICAgICAgIGlubmVyX2pvaW4oc3RhdGlvbnMsIGJ5ID0gYygiWm9uZSIsICJTdGF0aW9uIikpICU+JQ0KICAgICAgICBmaWx0ZXIoQ2VudGVycyA9PSBjZW50ZXIpICU+JQ0KICAgICAgICBnZ3Bsb3QoKSAlPiUNCiAgICAgICAgYXVzX2RhcmtfbWFwKCkgKyANCiAgICAgICAgZ2VvbV9wb2ludChhZXMoeCA9IExvbmdpdHVkZSwgeSA9IExhdGl0dWRlLCBjb2xvciA9IEdyb3VwKSkgKw0KICAgICAgICBzY2FsZV9jb2xvdXJfZGlzdGlsbGVyKHR5cGUgPSAiZGl2IiwgcGFsZXR0ZT0iUmRCdSIpICsNCiAgICAgICAgZ3VpZGVzKGNvbG91cj1GQUxTRSkNCiAgICBwcmludChwKQ0KfQ0KDQpgYGANCg0KIyMgUmVnaW9uYWwgR3JvdXBpbmcgQ29uY2x1c2lvbg0KDQpTRHMgYWxvbmUgaXQgbm90IHNvIGRpc3RpbmN0IGFuZCB0aGUgYWxnb3IgaGFzIHRyb3VibGUuIFVzZSBqdXN0IHRoZSBtZWFucyBzbyB0aGVpciBpcyBubyBpc3N1ZSBpbiBnZXR0aW5nIHZhbHVlcyBvbiB0aGUgc2FtZSBzY2FsZS4NCg0KRmluZHMgZ3JvdXAgd2l0aCBzaW1pbGFyIHJhaW5mYWxsIGRpc3RyaWJ1dGlvbnMuIFBpY2tzIG91dCB0aGUgYXJpZCByZWdpb25zLCBjb2FzdCByZWdpb25zIGFuZCBzbyBvbmUuDQoNCkhvd2V2ZXIgdGVuZHMgdG8gZG8gdGhpbmdzIGxpa2UgZ3JvdXAgYXJlYXMgb24gdGhlIGVhc3QgYW5kIHdlc3QgY29hc3Qgb2YgQXVzIHdoaWNoIEkga25vdyB0ZW5kIHRvIHZhcnkgc2VwZXJhdGVseS4NCg0KMTAgQ2VudGVycyBsb29rcyBhYm91dCBuaWNlLiBJdCBwaWNrcyBvdXQgdGhlIGFyZWFzIHRoYXQgYXJlIGluIGRyb3VnaHQgcmlnaHQgbm93Lg0KDQpTYXZlIGFzIHN0YXRpb25fcmVnaW9uYWxfY2x1c3RlcnMNCg0KYGBge3J9DQoNCnN0YXRpb25fcmVnaW9uYWxfY2x1c3RlcnMgPC0gc3RhdGlvbl9yYWluZmFsbF9tZWFucyAlPiUgDQogICAgc2VsZWN0KC1ab25lLCAtU3RhdGlvbikgJT4lIA0KICAgIGttZWFucygxMCkNCnN0YXRpb25fcmVnaW9ucyA8LSBzdGF0aW9uX3JhaW5mYWxsX21lYW5zDQpzdGF0aW9uX3JlZ2lvbnMkUmVnaW9uSUQgPC0gc3RhdGlvbl9yZWdpb25hbF9jbHVzdGVycyRjbHVzdGVyDQpzdGF0aW9uX3JlZ2lvbnMgPC0gc3RhdGlvbl9yZWdpb25zICU+JSBzZWxlY3QoWm9uZSwgU3RhdGlvbiwgUmVnaW9uSUQsIGV2ZXJ5dGhpbmcoKSkNCg0Kc3RhdGlvbnMgJT4lDQogICAgaW5uZXJfam9pbihzdGF0aW9uX3JlZ2lvbnMsIGJ5ID0gYygiWm9uZSIsICJTdGF0aW9uIikpICU+JQ0KICAgIGdncGxvdCgpICU+JQ0KICAgIGF1c19kYXJrX21hcCgpICsgDQogICAgZ2VvbV9wb2ludChhZXMoeCA9IExvbmdpdHVkZSwgeSA9IExhdGl0dWRlLCBjb2xvciA9IFJlZ2lvbklEKSkgICsNCiAgICBzY2FsZV9jb2xvdXJfZGlzdGlsbGVyKHR5cGUgPSAiZGl2IiwgcGFsZXR0ZT0iUmRCdSIpDQoNCg0KYGBgDQoNCiMjIFJlZ2lvbmFsIEdyb3VwIEFuYW1vbHkNCg0KR2VuZXJhdGUgYW5hbW9seSBncm91cGVkIGJ5IHJlZ2lvbnMNCg0KIyMjIENvbWJpbmVkIHBsb3QNCg0KYGBge3J9DQoNCnBsb3QodmFsdWF0aW9ucyAlPiUNCiAgICBpbm5lcl9qb2luKHN0YXRpb25zLCBieSA9IGMoIlpvbmUiLCAiU3RhdGlvbiIpKSAlPiUNCiAgICBpbm5lcl9qb2luKHN0YXRpb25fcmVnaW9ucywgYnkgPSBjKCJab25lIiwgIlN0YXRpb24iKSkgJT4lDQogICAgZmlsdGVyKE1ham9yX1BlcmlvZCA+IDE5MDApICU+JQ0KICAgIGdyb3VwX2J5KFJlZ2lvbklELCBNZWRpYW5fRGF0ZSkgJT4lDQogICAgc3VtbWFyaXNlKFJhaW5mYWxsLkFuYW1vbHkgPSBtZWFuKFJhaW5mYWxsLkFuYW1vbHkpKSAlPiUNCiAgICBnZ3Bsb3QoYWVzKHggPSBNZWRpYW5fRGF0ZSwgeSA9IFJhaW5mYWxsLkFuYW1vbHkpKSArDQogICAgZ2VvbV9saW5lKGNvbG9yPSJsaWdodGdyZXkiKSArDQogICAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxvZXNzIiwgZm9ybXVsYSA9IHkgfiB4LCBzcGFuID0gLjEpICsNCiAgICBmYWNldF93cmFwKH4gUmVnaW9uSUQpICsgeWxpbSgtMSwgMSkpDQoNCmBgYA0KDQojIyMgUGFnZWQgUGxvdA0KDQpgYGB7cn0NCg0KZm9yIChyZWdpb25JRCBpbiAxOjEwKSB7DQogICAgDQogICAgZ3JpZEV4dHJhOjpncmlkLmFycmFuZ2UoDQogICAgICAgIHN0YXRpb25zICU+JQ0KICAgICAgICAgICAgaW5uZXJfam9pbihzdGF0aW9uX3JlZ2lvbnMsIGJ5ID0gYygiWm9uZSIsICJTdGF0aW9uIikpICU+JQ0KICAgICAgICAgICAgZmlsdGVyKFJlZ2lvbklEID09IHJlZ2lvbklEKSAlPiUNCiAgICAgICAgICAgIGdncGxvdCgpICU+JQ0KICAgICAgICAgICAgYXVzX2RhcmtfbWFwKCkgKyANCiAgICAgICAgICAgIGdlb21fcG9pbnQoYWVzKHggPSBMb25naXR1ZGUsIHkgPSBMYXRpdHVkZSkpICsNCiAgICAgICAgICAgIGd1aWRlcyhjb2xvdXI9RkFMU0UpLA0KICAgICAgICANCiAgICAgICAgdmFsdWF0aW9ucyAlPiUNCiAgICAgICAgICAgIGlubmVyX2pvaW4oc3RhdGlvbnMsIGJ5ID0gYygiWm9uZSIsICJTdGF0aW9uIikpICU+JQ0KICAgICAgICAgICAgaW5uZXJfam9pbihzdGF0aW9uX3JlZ2lvbnMsIGJ5ID0gYygiWm9uZSIsICJTdGF0aW9uIikpICU+JQ0KICAgICAgICAgICAgZmlsdGVyKFJlZ2lvbklEID09IHJlZ2lvbklEKSAlPiUNCiAgICAgICAgICAgIGZpbHRlcihNYWpvcl9QZXJpb2QgPiAxOTAwKSAlPiUNCiAgICAgICAgICAgIGdyb3VwX2J5KE1lZGlhbl9EYXRlKSAlPiUNCiAgICAgICAgICAgIHN1bW1hcmlzZShSYWluZmFsbC5BbmFtb2x5ID0gbWVhbihSYWluZmFsbC5BbmFtb2x5KSkgJT4lDQogICAgICAgICAgICBnZ3Bsb3QoYWVzKHggPSBNZWRpYW5fRGF0ZSwgeSA9IFJhaW5mYWxsLkFuYW1vbHkpKSArDQogICAgICAgICAgICBnZW9tX2xpbmUoY29sb3I9ImdyZXkiKSArDQogICAgICAgICAgICBnZW9tX3Ntb290aChtZXRob2QgPSAibG9lc3MiLCBmb3JtdWxhID0geSB+IHgsIHNwYW4gPSAuMikgKw0KICAgICAgICAgICAgIHlsaW0oLTEsIDEpLA0KICAgICAgICANCiAgICAgICAgbmNvbCA9IDINCiAgICApDQp9DQoNCg0KYGBgDQoNCiMjIENsaW1hdGUgR3JvdXBpbmcgQ29uY2x1c2lvbg0KDQpJdCBzb3J0YSB3b3Jrcy4gWW91IGNhbiBzZWUgdGhlIG1ham9yIGRyb3VnaHQgaW4gZWFzdGVybiBhdXN0cmFsaWEuIEhvd2V2ZXIgc29tZSBhcmVhcyBhcmUgYXNzb2NpYXRlZCB0aGF0IHByb2JhYmx5IHNob3VsZG4ndCBiZS4NCg0KIyBUcnkgQ29ycmVsYXRpdmUgQ2x1c3RlcmluZw0KDQpDbHVzdGVyIHN0YXRpb25zIHRoYXQgdGVuZCB0byB2YXJ5IHRvZ2V0aGVyLg0KDQpUcnkga21lYW5zIHRvIHN0YXJ0IGNvcyBpdHMgZWFzeSENCg0KVHJ5IHNwcmVhZGluZyBhbGwgcGVyaW9kcyBkdXJpbmcgdGhlIGJhc2VsaW5lIHBlcmlvZCBpbnRvIG9uZSBzZXQgYW5kIHVzZSBrbWVhbnMgdG8gZ3JvdXAgc3RhdGlvbnMuDQoNCmBgYHtyfQ0Kc3RhdGlvbl9yYWluZmFsbF9hbmFtb2x5IDwtIHZhbHVhdGlvbnMgJT4lIA0KICAgIGZpbHRlcihNYWpvcl9QZXJpb2QgPj0gYXVzX2Jhc2VsaW5lX3N0YXJ0LCBNYWpvcl9QZXJpb2QgPD0gYXVzX2Jhc2VsaW5lX2VuZCkgJT4lDQogICAgbXV0YXRlKEJhc2VsaW5lX1BlcmlvZCA9IChNYWpvcl9QZXJpb2QgLSBhdXNfYmFzZWxpbmVfc3RhcnQpICogMTIgKyBNaW5vcl9QZXJpb2QpICU+JQ0KICAgIHNlbGVjdChab25lLCBTdGF0aW9uLCBBbmFtb2x5ID0gQmFzZWxpbmVfUGVyaW9kLCBSYWluZmFsbC5BbmFtb2x5KSAlPiUgDQogICAgc3ByZWFkKGtleSA9IEFuYW1vbHksIHZhbHVlID0gUmFpbmZhbGwuQW5hbW9seSwgc2VwID0gIi4iKSAlPiUNCiAgICB1bmdyb3VwKCkNCmBgYA0KDQojIyBHcm91cCBzdGF0aW9ucyBieSBhbmFtb2x5DQoNCmBgYHtyfQ0KDQpzdGF0aW9uX3JhaW5mYWxsX2FuYW1vbHlfY2x1c3RlcnMgPC0gZnVuY3Rpb24oY2VudGVycykgew0KICAgIHN0YXRpb25fa21lYW5zIDwtIHN0YXRpb25fcmFpbmZhbGxfYW5hbW9seSAlPiUgc2VsZWN0KC1ab25lLCAtU3RhdGlvbikgJT4lIGttZWFucyhjZW50ZXJzKQ0KICAgIHIgPC0gc3RhdGlvbl9yYWluZmFsbF9hbmFtb2x5DQogICAgciRDZW50ZXJzIDwtIGNlbnRlcnMNCiAgICByJFJlZ2lvbklEIDwtIHN0YXRpb25fa21lYW5zJGNsdXN0ZXINCiAgICByIDwtIHIgJT4lIHNlbGVjdChab25lLCBTdGF0aW9uLCBDZW50ZXJzLCBSZWdpb25JRCkNCiAgICByDQp9DQoNCmFuYW9tbHlfY2VudGVycyA8LSA0OjMwDQpzdGF0aW9uX3JhaW5mYWxsX2FuYW1vbHlfY2VudGVycyA8LSBhbmFvbWx5X2NlbnRlcnMgJT4lIG1hcF9kZnIoc3RhdGlvbl9yYWluZmFsbF9hbmFtb2x5X2NsdXN0ZXJzKQ0KDQpmb3IgKGNlbnRlciBpbiBhbmFvbWx5X2NlbnRlcnMpIHsNCiAgICBwbG90KHN0YXRpb25fcmFpbmZhbGxfYW5hbW9seV9jZW50ZXJzICU+JQ0KICAgICAgICBpbm5lcl9qb2luKHN0YXRpb25zLCBieSA9IGMoIlpvbmUiLCAiU3RhdGlvbiIpKSAlPiUNCiAgICAgICAgZmlsdGVyKENlbnRlcnMgPT0gY2VudGVyKSAlPiUNCiAgICAgICAgZ2dwbG90KCkgJT4lDQogICAgICAgIGF1c19kYXJrX21hcCgpICsgDQogICAgICAgIGdlb21fcG9pbnQoYWVzKHggPSBMb25naXR1ZGUsIHkgPSBMYXRpdHVkZSwgY29sb3IgPSBSZWdpb25JRCkpICsNCiAgICAgICAgc2NhbGVfY29sb3VyX2Rpc3RpbGxlcih0eXBlID0gImRpdiIsIHBhbGV0dGU9IlJkQnUiKSArDQogICAgICAgIGd1aWRlcyhjb2xvdXI9RkFMU0UpICsNCiAgICAgICAgZ2d0aXRsZShwYXN0ZTAoIkNlbnRlcnMgPSAiLCBjZW50ZXIpKSkNCn0NCg0KYGBgDQoNCiMjIyBQYWdlZCBQbG90DQoNCmBgYHtyfQ0KDQpzZWxlY3RlZF9hbmFvbWx5X2NlbnRlcnMgPC0gMTUNCnNlbGVjdGVkX3N0YXRpb25fcmFpbmZhbGxfYW5hbW9seSA8LSBzdGF0aW9uX3JhaW5mYWxsX2FuYW1vbHlfY2VudGVycyAlPiUgZmlsdGVyKENlbnRlcnMgPT0gc2VsZWN0ZWRfYW5hb21seV9jZW50ZXJzKQ0Kc2VsZWN0ZWRfYW5hb21seV9zdGFydCA8LSAxOTUwDQpzZWxlY3RlZF9hbmFvbWx5X2VuZCA8LSAyMDE4DQoNCmZvciAocmVnaW9uSUQgaW4gMTpzZWxlY3RlZF9hbmFvbWx5X2NlbnRlcnMpIHsNCiAgICANCiAgICBncmlkRXh0cmE6OmdyaWQuYXJyYW5nZSgNCiAgICAgICAgc3RhdGlvbnMgJT4lDQogICAgICAgICAgICBpbm5lcl9qb2luKHNlbGVjdGVkX3N0YXRpb25fcmFpbmZhbGxfYW5hbW9seSwgYnkgPSBjKCJab25lIiwgIlN0YXRpb24iKSkgJT4lDQogICAgICAgICAgICBmaWx0ZXIoUmVnaW9uSUQgPT0gcmVnaW9uSUQpICU+JQ0KICAgICAgICAgICAgZ2dwbG90KCkgJT4lDQogICAgICAgICAgICBhdXNfZGFya19tYXAoKSArIA0KICAgICAgICAgICAgZ2VvbV9wb2ludChhZXMoeCA9IExvbmdpdHVkZSwgeSA9IExhdGl0dWRlKSkgKw0KICAgICAgICAgICAgZ3VpZGVzKGNvbG91cj1GQUxTRSksDQogICAgICAgIA0KICAgICAgICB2YWx1YXRpb25zICU+JQ0KICAgICAgICAgICAgaW5uZXJfam9pbihzdGF0aW9ucywgYnkgPSBjKCJab25lIiwgIlN0YXRpb24iKSkgJT4lDQogICAgICAgICAgICBpbm5lcl9qb2luKHNlbGVjdGVkX3N0YXRpb25fcmFpbmZhbGxfYW5hbW9seSwgYnkgPSBjKCJab25lIiwgIlN0YXRpb24iKSkgJT4lDQogICAgICAgICAgICBmaWx0ZXIoUmVnaW9uSUQgPT0gcmVnaW9uSUQpICU+JQ0KICAgICAgICAgICAgZmlsdGVyKE1ham9yX1BlcmlvZCA+PSBzZWxlY3RlZF9hbmFvbWx5X3N0YXJ0LCBNYWpvcl9QZXJpb2QgPD0gc2VsZWN0ZWRfYW5hb21seV9lbmQpICU+JQ0KICAgICAgICAgICAgZ3JvdXBfYnkoTWVkaWFuX0RhdGUpICU+JQ0KICAgICAgICAgICAgc3VtbWFyaXNlKFJhaW5mYWxsLkFuYW1vbHkgPSBtZWFuKFJhaW5mYWxsLkFuYW1vbHkpKSAlPiUNCiAgICAgICAgICAgIGdncGxvdChhZXMoeCA9IE1lZGlhbl9EYXRlLCB5ID0gUmFpbmZhbGwuQW5hbW9seSkpICsNCiAgICAgICAgICAgIGdlb21fbGluZShjb2xvcj0iZ3JleSIpICsNCiAgICAgICAgICAgIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsb2VzcyIsIGZvcm11bGEgPSB5IH4geCwgc3BhbiA9IC4xKSwNCiAgICAgICAgDQogICAgICAgIG5jb2wgPSAyDQogICAgKQ0KfQ0KDQoNCmBgYA0K